home *** CD-ROM | disk | FTP | other *** search
/ Collection of Tools & Utilities / Collection of Tools and Utilities.iso / pascal / tpl60n19.zip / TESTPRGS.ZIP / DEXP.PAS < prev    next >
Pascal/Delphi Source File  |  1993-02-14  |  8KB  |  292 lines
  1. PROGRAM DExp;    { ported from Fortran original 05-01-92 Norbert Juffa }
  2.  
  3. {$A+,B-,D-,E+,F-,G-,I-,L-,N-,O-,R-,S-,V-,X-}
  4.  
  5. USES MachArit, Power;
  6.  
  7. {     PROGRAM TO TEST DEXP
  8. C
  9. C     DATA REQUIRED
  10. C
  11. C        NONE
  12. C
  13. C     SUBPROGRAMS REQUIRED FROM THIS PACKAGE
  14. C
  15. C        MACHAR - AN ENVIRONMENTAL INQUIRY PROGRAM PROVIDING
  16. C                 INFORMATION ON THE FLOATING-POINT ARITHMETIC
  17. C                 SYSTEM.  NOTE THAT THE CALL TO MACHAR CAN
  18. C                 BE DELETED PROVIDED THE FOLLOWING FOUR
  19. C                 PARAMETERS ARE ASSIGNED THE VALUES INDICATED
  20. C
  21. C                 IBETA - THE RADIX OF THE FLOATING-POINT SYSTEM
  22. C                 IT    - THE NUMBER OF BASE-IBETA DIGITS IN THE
  23. C                         SIGNIFICAND OF A FLOATING-POINT NUMBER
  24. C                 XMIN  - THE SMALLEST NON-VANISHING FLOATING-POINT
  25. C                         POWER OF THE RADIX
  26. C                 XMAX  - THE LARGEST FINITE FLOATING-POINT NO.
  27. C
  28. C        REN(K) - A FUNCTION SUBPROGRAM RETURNING RANDOM REAL
  29. C                 NUMBERS UNIFORMLY DISTRIBUTED OVER (0,1)
  30. C
  31. C
  32. C     STANDARD FORTRAN SUBPROGRAMS REQUIRED
  33. C
  34. C         DABS, DINT, DLOG, DMAX1, DEXP, DFLOAT, DSQRT
  35. C
  36. C
  37. C     LATEST REVISION - DECEMBER 6, 1979
  38. C
  39. C     AUTHOR - W. J. CODY
  40. C              ARGONNE NATIONAL LABORATORY
  41. C
  42. }
  43.  
  44.  
  45.  
  46. FUNCTION REN (K: LONGINT): REAL;
  47.  
  48. {
  49.       DOUBLE PRECISION FUNCTION REN(K)
  50. C
  51. C     RANDOM NUMBER GENERATOR - BASED ON ALGORITHM 266 BY PIKE AND
  52. C      HILL (MODIFIED BY HANSSON), COMMUNICATIONS OF THE ACM,
  53. C      VOL. 8, NO. 10, OCTOBER 1965.
  54. C
  55. C     THIS SUBPROGRAM IS INTENDED FOR USE ON COMPUTERS WITH
  56. C      FIXED POINT WORDLENGTH OF AT LEAST 29 BITS.  IT IS
  57. C      BEST IF THE FLOATING POINT SIGNIFICAND HAS AT MOST
  58. C      29 BITS.
  59. C
  60. }
  61.  
  62. VAR   J:  LONGINT;
  63. CONST IY: LONGINT = 100001;
  64.  
  65. BEGIN
  66.    J  := K;
  67.    IY := IY * 125;
  68.    IY := IY - (IY DIV 2796203) * 2796203;
  69.    REN:= 1.0 * (IY) / 2796203.0e0 * (1.0e0 + 1.0e-6 + 1.0e-12);
  70. END;
  71.  
  72.  
  73.  
  74. FUNCTION MAX1 (A, B:REAL): REAL;
  75. BEGIN
  76.    IF A > B THEN
  77.       MAX1 := A
  78.    ELSE
  79.       MAX1 := B;
  80. END;
  81.  
  82.  
  83.  
  84. VAR   I,IBETA,IEXP,IOUT,IRND,IT,I1,J,K1,K2,K3,MACHEP,
  85.       MAXEXP,MINEXP,N,NEGEP,NGRD: LONGINT;
  86.       A,AIT,ALBETA,B,BETA,D,DEL,EPS,EPSNEG,ONE,R6,R7,
  87.       TWO,TEN,V,W,X,XL,XMAX,XMIN,XN,X1,Y,Z,ZERO,ZZ,
  88.       HALF,NINETENTH,FOUR,FORTYFIVE,SIXTEEN,SIXTEENTH,
  89.       SIXHUNDREDTH: REAL;
  90.  
  91. LABEL 100, 110, 120, 270, 300;
  92.  
  93.  
  94. BEGIN
  95.  
  96.    N := 1000000;    { number of trials }
  97.  
  98.    MACHAR (IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP,MAXEXP,
  99.            EPS,EPSNEG,XMIN,XMAX);
  100.    PRINTPARAM (IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP,MAXEXP,
  101.                EPS,EPSNEG,XMIN,XMAX);
  102.    BETA         := IBETA;
  103.    ALBETA       := LN (BETA);
  104.    AIT          := IT;
  105.    ZERO         := 0;
  106.    ONE          := 1;
  107.    TWO          := 2;
  108.    FOUR         := 4;
  109.    TEN          := 10;
  110.    SIXTEEN      := 16;
  111.    FORTYFIVE    := 45;
  112.    HALF         := 0.5;
  113.    NINETENTH    := 0.9;
  114.    SIXTEENTH    := 0.0625;
  115.    SIXHUNDREDTH := 0.006;
  116.    V            := SIXTEENTH;
  117.    A            := TWO;
  118.    B            := LN (A) * HALF;
  119.    A            := -B + V;
  120.    D            := LN (NINETENTH*XMAX);
  121.    XN           := N;
  122.    I1           := 0;
  123.  
  124. {---------------------------------------------------------------------}
  125. {     RANDOM ARGUMENT ACCURACY TESTS                                  }
  126. {---------------------------------------------------------------------}
  127.  
  128.    FOR J := 1 TO 3 DO BEGIN
  129.       K1 := 0;
  130.       K3 := 0;
  131.       X1 := ZERO;
  132.       R6 := ZERO;
  133.       R7 := ZERO;
  134.       DEL:= (B - A) / XN;
  135.       XL := A;
  136.  
  137.       FOR I := 1 TO N DO BEGIN
  138.          X := DEL * REN(I1) + XL;
  139.  
  140. {---------------------------------------------------------------------}
  141. {     PURIFY ARGUMENTS                                                }
  142. {---------------------------------------------------------------------}
  143.  
  144.          Y := X - V;
  145.          IF Y < ZERO THEN
  146.             X := Y + V;
  147.          Z  := EXP (X);
  148.          ZZ := EXP (Y);
  149.          IF J = 1 THEN
  150.             GOTO 100;
  151.          IF IBETA <> 10 THEN
  152.             Z := Z * SIXTEENTH - Z *
  153.                  2.4453321046920570389e-3; { 1/16 - exp (-45/16) }
  154.          IF IBETA = 10 THEN
  155.             Z := Z * SIXHUNDREDTH +  Z *
  156.                  5.466789530794296106e-5;  { 6/100 - exp (-45/16) }
  157.          GOTO 110;
  158. 100:     Z := Z - Z * 6.058693718652421388e-2;  { 1 - exp (-1/16) }
  159. 110:     IF Z <> ZERO THEN
  160.             W := (Z - ZZ) / Z
  161.          ELSE IF ZZ <> 0 THEN
  162.             W := ONE;
  163.          IF W > ZERO THEN
  164.             K1 := K1 + 1;
  165.          IF W < ZERO THEN
  166.             K3 := K3 + 1;
  167.          W := ABS (W);
  168.          IF W <= R6 THEN
  169.             GOTO 120;
  170.          R6 := W;
  171.          X1 := X;
  172. 120:     R7 := R7 + W * W;
  173.          XL := XL + DEL;
  174.       END;
  175.  
  176.       K2 := N - K3 - K1;
  177.       R7 := SQRT (R7/XN);
  178.  
  179.       WRITELN;
  180.       WRITELN;
  181.       WRITELN ('TEST OF EXP (X-', V:15, ') VS EXP(X)/EXP(', V:15, ')');
  182.       WRITELN;
  183.       WRITELN (N, ' RANDOM ARGUMENTS WERE TESTED FROM THE INTERVAL');
  184.       WRITELN ('(', A, ',', B, ')');
  185.       WRITELN;
  186.       WRITELN ('EXP(X-V) WAS LARGER', K1:6, ' TIMES');
  187.       WRITELN ('             AGREED', K2:6, ' TIMES');
  188.       WRITELN ('    AND WAS SMALLER', K3:6, ' TIMES');
  189.       WRITELN;
  190.       WRITELN ('THERE ARE ', IT, ' BASE ', IBETA,
  191.                ' SIGNIFICANT DIGITS IN A FLOATING-POINT NUMBER');
  192.       WRITELN;
  193.       W := -999;
  194.       IF R6 <> ZERO THEN
  195.          W := LN (ABS(R6))/ALBETA;
  196.       WRITELN ('THE MAXIMUM RELATIVE ERROR OF          ', R6:12,
  197.                ' = ', IBETA, ' **', W:7:2);
  198.       WRITELN ('OCCURED FOR X = ', X1);
  199.       W := MAX1 (AIT+W,ZERO);
  200.       WRITELN;
  201.       WRITELN ('THE ESTIMATED LOSS OF BASE ', IBETA,
  202.                ' SIGNIFICANT DIGITS IS        ', W:7:2);
  203.       W := -999.0;
  204.       IF R7 <> ZERO THEN
  205.          W := LN (ABS(R7))/ALBETA;
  206.       WRITELN;
  207.       WRITELN ('THE ROOT MEAN SQUARE RELATIVE ERROR WAS', R7:12,
  208.                ' = ', IBETA, ' **' , W:7:2);
  209.       W := MAX1 (AIT+W,ZERO);
  210.       WRITELN;
  211.       WRITELN ('THE ESTIMATED LOSS OF BASE ', IBETA,
  212.                ' SIGNIFICANT DIGITS IS        ', W:7:2);
  213.       IF J = 2 THEN
  214.          GOTO 270;
  215.       V := FORTYFIVE / SIXTEEN;
  216.       A := -TEN * B;
  217.       B := FOUR * XMIN * POW (BETA ,IT);
  218.       B := LN (B);
  219.       GOTO 300;
  220. 270:  A := -TWO * A;
  221.       B := TEN * A;
  222.       IF B < D THEN
  223.          B := D;
  224. 300:
  225.    END;
  226.  
  227. {---------------------------------------------------------------------}
  228. {    SPECIAL TESTS                                                    }
  229. {---------------------------------------------------------------------}
  230.  
  231.    WRITELN;
  232.    WRITELN;
  233.    WRITELN ('SPECIAL TESTS');
  234.    WRITELN;
  235.    WRITELN ('THE IDENTITY EXP(X)*EXP(-X) = 1.0 WILL BE TESTED');
  236.    WRITELN;
  237.    WRITELN ('          X           F(X)*F(-X)-1');
  238.  
  239.    FOR I := 1 TO 5 DO BEGIN
  240.       X := REN(I1) * BETA;
  241.       Y := -X;
  242.       Z := EXP(X) * EXP(Y);
  243.       Z := Z - ONE;
  244.       WRITELN (X:18, Z:18);
  245.    END;
  246.  
  247.    WRITELN;
  248.    WRITELN;
  249.    WRITELN ('TEST OF SPECIAL ARGUMENTS');
  250.    X := ZERO;
  251.    Y := EXP (X) - ONE;
  252.    WRITELN;
  253.    WRITELN ('EXP (0.0) - 1.0 = ', Y:15);
  254.    X := INT (LN(XMIN));
  255.    Y := EXP (X);
  256.    WRITELN;
  257.    WRITELN ('EXP (', X:15, ') = ', Y:15);
  258.    X := INT (LN(XMAX));
  259.    Y := EXP (X);
  260.    WRITELN;
  261.    WRITELN ('EXP (', X:15, ') = ', Y:15);
  262.    X := X / TWO;
  263.    V := X / TWO;
  264.    Y := EXP (X);
  265.    Z := EXP (V);
  266.    Z := Z * Z;
  267.    WRITELN;
  268.    WRITELN ('IF EXP (', X:15, ') = ', Y:15, ' IS NOT ABOUT ');
  269.    WRITELN ('EXP (', V:15, ')**2 = ', Z:15, ' THERE IS AN ARG RED ERROR');
  270.  
  271. {---------------------------------------------------------------------}
  272. {     TEST OF ERROR RETURNS                                           }
  273. {---------------------------------------------------------------------}
  274.  
  275.    WRITELN;
  276.    WRITELN;
  277.    WRITELN ('TEST OF ERROR RETURNS');
  278.    WRITELN;
  279.    X := -ONE / SQRT (XMIN);
  280.    WRITELN ('EXP WILL BE CALLED WITH THE ARGUMENT ', X:15);
  281.    Y := EXP (X);
  282.    WRITELN ('EXP RETURNED THE VALUE', Y:15);
  283.    WRITELN;
  284.    X := -X;